Viscous Flow and Jump Dynamics in Molecular Supercooled Liquids: 
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The transport and the relaxation properties of a molecular supercooled liquid on a isobar is studied 
by molecular dynamics. The molecule is a rigid heteronuclear biatomic system. The diffusivity 
is fitted over four orders of magnitude by the power law D oc (T — T C ) 1D with yr> = 1.93 ± 0.02 
and T c = 0.458 ± 0.002. The self-part of the intermediate scattering function F s (k rnax ,t) exhibits 
a step-like behavior at the lowest temperatures. On cooling, the increase of the related relaxation 
time t q tracks the diffusivity, i.e. T a oc (k^^D) -1 . At the lowest temperatures fractions of highly 
mobile and trapped molecules are evidenced. Translational jumps are evidenced. The duration of the 
jumps exhibits a distribution. The distribution of the waiting-times before a jump takes place ip(t) 
is exponential at higher temperatures. At lower temperatures a power-law divergence is evidenced 
at short times, %jj(t) oc t^^ 1 with < £ < 1 which is ascribed to intermittency. The shear viscosity 
is fitted by the power law r\ oc (T — T c ) lri with 7^ = —2.20 ± 0.03 at the lowest temperatures. At 
higher temperatures the Stokes-Einstein fits the data if stick boundary conditions are assumed. The 
product Dn/T increases at lower temperatures and the Stokes-Einstein relation breaks down at a 
temperature which is close to the one where the intermittency is evidenced by tp(t). A precursor 
effect of the breakdown is observed which manifests as an apparent stick-slip transition. 

PACS numbers: 64.70.Pf, 02.70.Ns, 66.20.+d , 66.10.-x 



I. INTRODUCTION 

The relaxation phenomena and the tranport properties 
of supercooled liquids and glassy materials are topics of 
current interest p],0]. It is well known that, on approach- 
ing the glass transition temperature T g from above, diffu- 
sion coefficients and relaxation times exhibit remarkable 
changes of several orders of magnitude which are under 
intense experimental, theoretical and numerical investi- 
gation. In the high-temperature regime the changes usu- 
ally track the shear viscosity r\ in the sense that , if X 
denotes the diffusivity or the inverse of a relaxation time, 
the product Xr]/T is nearly temperature- independent. 
In particular, both the Stokes-Einstein, D w kT/6nria, 
and the Debye-Stokes-Einstein laws, D r sa kT/r/a 3 are 
found to work nicely, D, D r and a being the transla- 
tional and the rotational diffusivity and the molecular 
radius, respectively. Differently, in deeply supercooled 
regimes there is wide evidence that the product increases 
on cooling evidencing the breakdown of the hydrodinamic 
behavior at molecular level and the decoupling by the 
viscous flow HH [j}{iJl. 

Several models suggest that the decoupling between 
microscopic time scales and the viscous flow is a sig- 
nature of the heterogeneous dynamics close to the glass 
transition, i.e. a spatial distribution of transport and re- 
laxation properties [^6[-p9|. Interesting alternatives are 
provided by frustrated lattice gas models j2(| and the 
"energy landscape" picture |^21-24|. Most intepreta- 
tions suggest the existence of crossover temperatures be- 
low which a change of relaxation mechanism must occur 



[pj. These temperatures are broadly found around 1.2T g , 
i.e. in the region where the critical temperature T c pre- 
dicted by the mode-coupling theory of the glass transition 



( MCT ) is found |25). Recent extensions of MCT for the 
shear viscosity to take into account current-fluctuations 
in addition to density fluctuations are reported |2(J . 

During the last years molecular dynamics simulations 
( MD ) proved to be a powerful tool to investigate su- 
percooled liquids ( for a recent review see ref. [^?J ). 
To date, MD studies investigated decoupling phenomena 
in atomic pure liquids [ p4j and atomic binary mixtures 
p8| [33[ . Most MD studies confirmed that the decoupling 
is due to dynamic heterogeneities [ p8pC| |34]|. In fact, 
"active" [BlJ or "mobile" pi] regions which largely con- 
tribute to set the macroscopic average value have been 
identified. In such regions hopping processes, enhanc- 
ing the transport with respect to the hydrodynamic be- 
havior, have been evidenced p8| , pT[ . The occurrence of 
jumps in glasses has been reported several times in the 
recent past 

We are not aware of MD studies of the decoupling, and, 
more generally, of dynamic heterogeneities in molecular 
systems. This motivated the present and the following 
paper |38| to investigate the transport and the relax- 
ation in a three-dimensional, one-component, molecular 
system. This is an important feature since most exper- 
imental work is carried out on that class of materials. 
In particular, the decoupling of microscopic relaxation 
from the viscous flow will be addressed in the light of 
the increased role played by the hopping processes. The 
present paper is limited to the translational degrees of 
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freedom whereas the following paper [B8| deals with the 
rotational degrees of freedom. 

The paper is organized as follows. In Sec. ||the model 
and the details of the simulation are presented. The re- 
sults concerning both singl e-p article and collective prop- 
erties are discussed in Sec. 
summarized in Sec. 
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III . The main conclusions are 



6, 5, 3, 2, 1.4, 1.1, 0.85, 0.77, 0.70, 0.632, 0.588, 0.549, 0.52, 0.5. 
The characterization of the system at T = 0.77 was lim- 
ited to the diffusivity and the viscosity. The density at 
T = 0.5 was 0.6998 and decreased of a factor of about 
three at the highest temperatures. To check that thermal 
hystory is negligible, at least two independent equilibra- 
tions were performed and the subsequent production runs 
compared at the lowest temperatures. 



II. MODEL AND DETAILS OF SIMULATION 

The system under study is a biatomic molecular liquid. 
The model has been extensively investigated to test the 
MCT predictions The atoms A and B of each 

molecule have mass m and are spaced by d. Atoms on 
different molecules interact via the Lennard- Jones poten- 
tial: 

V a p(r) = 4e Q/3 (cr Q/3 /r) 12 - {a a p/r) 6 , a, f3 E {A,B} 

(1) 

The potential was cutoff and shifted at r cuto ff — 
2.5(taa- Henceforth, reduced units will be used. Lenghts 
are in units of gaa , energies in units of caa and masses in 

(ma 2 \ 1/2 

units of m. The time unit is ( J , corresponding 

to about 2ps for the Argon atom. The pressure P, tem- 
perature T and shear viscosity r\ are in units of caa/ ®aa 
, e AA /kB and ^fmeAA j 'cr AA , respectively. 

The model parameters in reduced units are: o aa — 
(TAB = 1.0, o B B = 0.95, 6AA = £ab = 1-0, e B B = 0.95, 
d = 0.5, rriA = tub = m = 1.0. The <jaa and ubb 
values were chosen to avoid crystallization. The sample 
has N = Nat/ '2 = 1000 molecules which are accommo- 
dated in a cubic box with periodic boundary conditions. 
The viscosity was evaluated by using samples of N = 108 
molecules. 

We examined the isobar at P = 1.5 by the follow- 
ing procedure. First, the sample was equilibrated in 
isothermal-isobaric conditions, for a time t eq t eq was at 
least one order of magnitude longer the time needed by 
the self-part of the intermediate scattering function eval- 
uated at the maximum of the static structure factor to 
become smaller than 0.1. In this step the equations of 
motion were integrated by using the RATTLE algorithm 
with Nose- Andersen constant temperature and pressure 
dynamics [Q. The algorithm is detailed in Appendix 
|A|. The starting conditions of the equilibration run make 
the total momentum of the system to vanish and locate 
the center of mass at the center of the box. The Nose- 
Andersen Lagrangian ensures that the center-of-mass po- 
sition and the total momentum will not change during the 
run. 

The data were collected in a production run in micro- 
canonical conditions. Integration was carried out by the 
RATTLE algorithm. The St step ranges from 0.001 at 
higher temperatures to 0.004 at lower ones in the pro- 
duction run. The temperatures we investigated are T = 



III. RESULTS AND DISCUSSION 

The section discusses the results of the study. First, we 
characterize the static properties of the system, i.e. the 
radial distribution function of the center-of-mass g(r) and 
the static structure factor S(k). Then, the single-particle 
and the collective dynamical properties of the system will 
be presented. 



A. Static properties 

The radial distribution function of the center of mass 
g(r) is defined as: 

5(r)= 47A^5: (,5(|R ' (o) - R ^ o) - r|)> (2) 

p is the average density and (t) is the position of the 
center of mass of the i— th molecule at time t. Represen- 
tative plots of g(r) at different temperatures are shown in 
Fig. . The pattern is typical of a disordered system. At 
T = 0.5, g(r) exhibits the maximum at r = 1.2 and the 
minimum at r = 1.6. On increasing the temperature, the 
peaks broaden and shift at higher r values. The shoulder 
which is observed at r = 1.1 at the lowest temperatures 
has been already observed E3| and may be ascribed to 
local "T-shaped" or cross configurations of the molecules 

0. 

The static structure factor S(k) is defined as 
1 



S(k) 



N 



(pkP-k 



(3) 



Pk is the Fourier transform of the density. Representative 
plots of S(k) at different temperatures are shown in fig. 
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FIG. 1. Radial distribution function g(r) ( top ) and static 
structure factor S(k) ( bottom ) of the molecular center of 
mass for selected temperatures. 

The above results rule out possible crystallizations of 
the sample. This is anticipated by the absence of anoma- 
lies in the temperature dependence of the density. Direct 
inspection of several snapshots of the molecular configu- 
rations also supported the conclusion and did not reveal 
any orientational ( liquid-crystalline ) order. The absence 
of orientational order was also evidenced by the long-time 
behaviour of the rotational correlation function |3§| ]. 

B. Single-particle dynamics 

The section will discuss the single-particle dynamics. 
At short time scales one quantity of interest is the veloc- 
ity self-correlation function of the center of mass : 



0»(t) = -jy v,(0)) (4) 

i=l 

(A(t)) denotes a proper time average of A(t). Fig. |^ 
plots C vv for different temperatures. 




t 

FIG. 2. Velocity self-correlation function of the center of 
mass. The curves refer to T = 3, 2, 1.4, 1.1, 0.85, 0.7, 0.5. 

At very short times C vv (t) = < v 2 > (1- < Q 2 > 
t 2 /2 + . . .), where < il 2 > is the Einstein frequency Q. 
At the lower temperatures the damped oscillatory motion 
of the molecule inside the cage where it is accommodated 
becomes apparent. The effect signals that on a short time 
scale the supercooled liquid behaves as a solid. 

The mean squared displacement of center of mass, R(t) 
provides a first view of the translational motion of the 
molecules at intermediate and long time scales. It is de- 
fined as: 

fl(*) = ^E<|Ri(*)-Ri(o)| 2 >. (5) 

1=1 

Plots of R(t) at different temperatures are shown in ng||. 
At short times the motion is ballistic and R(t) oc t 2 . At 
intermediate times R(t) exhibits a crossover regime which 
is strongly temperature dependent. At higher tempera- 
ture it reduces to a knee joining the ballistic and the 
diffusive regimes. At lower temperatures the molecule is 
effectively trapped inside the cage of the first neighbours 
and R(t) exhibits a plateau. At long times the motion is 
diffusive and R(t) = 6Dt where D is the diffusion con- 
stant. 
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FIG. 3. Mean squared displacement of molecular center-of- 
mass R(t) for all the temperatures under investigation but 
T = 0.77. 

The occurrence of oscillatory motion during trapping 
is shown explicitely in fig.^ where the mean squared dis- 
placement and the velocity correlation function are com- 
pared at T = 0.5. The plot evidences that velocity cor- 
relations got lost within the lifetime of the cage. 
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FIG. 4. Comparison of the mean squared displacement 
R(t) and the velocity correlation function C vv (t) at T = 0.5. 
The onset of the damped oscillations takes place at the be- 
ginning of the plateau of R(t): the molecule is colliding with 
the cage of the first neighbours. 

The translational diffusion coefficient D is evaluated 
by the Einstein relation (HQ : 

R(t) 



D = lim 
t- 



(it 



(6) 



In Fig. ^| the temperature dependence of D is shown. 
The plot emphasizes the scaling property of D which is 
nicely described over four decades by the power law : 



Theoretical justification of eq.0 is provided by MCT 
p5| . In particular, the ideal MCT predicts the inequality 
7d > 1.5. Our best fit values are = 1-93 ± 0.02 and 
T c = 0.458 ± 0.002, C D = 0.0481 ± 0.0004. For the same 
biatomic system with N = 500 and P = 1 it was found 
jn = 2.2 and T c = 0.475 M. 



1 2 3 4 5 
I i i i i I i i i i I i i i— i— I— i— i— i— i— I— i 



10° ^ 



10" - 



° IO" 2 -4 



10"" - 



T c =0.458 ±0.002 



-i — i — i 1 1 1 1 1 — 

3 4 5 6 1 

0.1 



- 1 — I — I I I 1 1 1 

3 4 5 6 1 



1 



- 1 — i — i r 

3 4 5 6 



T-T c 

FIG. 5. Temperature dependence of the translational dif- 
fusion coefficient D. The dashed line is a fit with the power 
law eqj| with -y D = 1.93 ± 0.02 , T c = 0.458 ± 0.002 and 
C D = 0.0481 ± 0.0004. 

It must be noted that at lower temperatures the above 
power law is expected to underestimate the diffusion co- 
efficient El]]. The larger diffusivity with respect to the 
ideal MCT prediction close to T c is not surprising. Ac- 
cording to this approach, density-density correlations do 
not vanish at long times below T c leading to a divergence 
of D at T c . The extended MCT points out that hop- 
ping processes provide an effective mechanism to relax 
the density fluctuations at long time and make the dif- 
fusivity finite close to T c . This point will be explicitely 
discussed later. 

To characterize further the single-particle dynamics 
we evaluated the self-part of the intermediate scattering 
function : 



F s (k,t) = (ex P [ik ■ (Rjit) - Rj(0))]) 



(8) 



In fig. H F s (k,t) is plotted for all the temperatures we 
investigated at k = k max . k max is the position of the 
main peak of the static structure factor S(k) at temper- 
ature T. We note that F s (k, t) always decays to zero and 
is independent of the thermal hystory of the sample ( i.e. 
independent runs lead to the same result ) signaling the 
equilibration of the sample. At higher temperatures the 
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decay is virtually exponential at long times, whereas at 
lower temperatures a two-step decay is observed. 




FIG. 6. Self-part of the intermediate scattering function 
F 3 (kmax,t). The curves refer to all the temperatures under 
investigation but T = 0.77. 

The plateau which is observed at lower temperatures 
in F s (k max ,t) is due to the trapping of the molecules. 
This is seen by comparing hg.|| with fig|| Precise predic- 
tions on the scaling features of the plateau of F s (k max ,t) 
are made by MCT p5| . A thorough comparison for the 
present model is found in ref. fjj]| . 
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FIG. 7. Temperature dependence of the product T Q _Dfc^ a:r . 

At longer times the cage where the molecule is trapped 
opens and the density-density correlations vanish. The 
long-time decay of F s (k max ,t) is well fitted by the 
stretched exponential Aexp(— t/r) 13 . Stretching is appre- 
ciable at lower temperatures ( /3o.5 = 0.677 ) and becomes 
negligible ( /3 w 1 ) for T > 0.85. 

The knowledge of F s (k max ,t) offers the opportunity 
to investigate the time scale of structural relaxation be- 
ing usually denoted by the a- relaxation time r a . We 



have evaluated it by the condition F s (k max , r Q ) = 1/e = 
0.3679. At low temperatures a nice fit is provided by 
i~ a = C a {T — T c ) 7 =, being T c the same best- fit value 
drawn by the diffusivity. The best fit value of the ex- 
ponent is 7 Q = —1.89 ± 0.05 which is quite close to 
7d = 1.93. Diffusive motion in simple liquids leads 
to the exponential decay of F s (k,t) with time costant 
t = l/Dk 2 [Q. Even the decay at long times is better 
described by a stretched exponential, the approximate 
equality j a = jjj prompted us to investigate the temper- 
ature dependence of the product T a Dk 2 nax . The results 
are shown in fig.0. For T < 1.1 the product approaches 
the constant value 0.72 ± 0.05. In this range r Q and D 
change of more than two orders of magnitude. Alter- 
native definitions of r a , e.g. involving the area below 
F s {k maxi t) do not affect significantly the result. Remov- 
ing the fc^ a2 . term results in sligthly poorer fit. The k 2 
scaling of the primary relaxation time t^, evaluated as 
the area below F s (k,t), was already noted in the hard 
sphere system at the critical packing fraction in the re- 
gion 1.5 < ka < 30 45 . In particular, it was found 
r^Dk 2 ~ 1.08 in good agreement with our result and 
F s (k,t) exhibits stretching with f3(k max ) ~ 0.8 . 

The above findings agree with some predictions of 
MCT: i) T c must be independent of the quantity under 
study, ii) j a must be equal to jd- The equality 7 Q = jz> 
must be taken with great caution. Other studies on the 
same model system with N — 500 and P = 1 found 
la > Id even if fitting the diffusivity and r Q provide the 
same T c value [Ell. 
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FIG. 8. Comparison of the non-gaussian parameter cti and 
the mean squared displacement R(t) at T — 0.5. The maxi- 
mum of the non-gaussian parameter occurs at t* . t* provides 
an estimate of the trapping time and is comparable with T a - 

Additional information on the single-particle dynamics 
is provided by the self part of the Van Hove function 
G.(r,t)@: 



1 N 

Gs(r,t) = -C£S(* + K t (0) - ^(t))) 



(9) 
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In isotropic liquids the Van Hove function depends only 
on the modulus r of r. The interpretation of G s (r,t) is 
direct. The product G s (r, t) ■ Anr 2 is the probability that 
the molecule is at a distance between r and r + dr from 
the initial position after a time t. 




FIG. 9. The ratio N 3 (r) = (G s (r,f) - Gf(r,t*)) 
/G 3 s (r,t*) plotted at T — 0.5(dashed line), T = 0.588 (solid 
line), T = 0.85 (dot-dashed line). Inset: same quantity on a 
logarithm scale. 

The shape of G s (r, t) in the system under study does 
not reveal particular features |41j| . In particular, differ- 
ently from other studies |^9| , no secondary peak at r « 1 , 
the nearest-neighbor distance, is found to support the 
conclusion that hopping is occurring. A deeper insight is 
gained by comparing G s (r, t) with the gaussian approxi- 
mation: 

G»(r, t) = [3/27r(r 2 (i))] 3 / 2 cxp(-3r 2 /2(r 2 (t)» (10) 

Eq. is the correct limit of G s (r, t) at short ( ballistic 
regime, (r 2 (t)) = 3kT/mt 2 ) and long times ( diffusion 
regime, (r 2 (t)) = 6Dt ) [j44f . Thirumalai and Mountain 
p8[ and more recently Kob et al. Q noted discrepancies 
between G s (r, t) and G§ (r, t) in supercooled atomic mix- 
tures. It is found that G s (r,t) exceeds G a a {r,t) at short 
and long r values. The effect is particularly evident at 
t w t* where t* is the time where the non-gaussian pa- 
rameter 



a a (t) = 3(r\t))/5(r 2 (t)) 2 - 1, 



(11) 



reaches the maximum value EJ. At t* the stochastic 
properties of r differ by the ones of a gaussian variable 
to the maximum. A plot of oi2(t) is shown in fig|| where 
it is compared to the mean squared displacement R{t) 
at T — 0.5. The maximum of ot2(t) at t = t* is located 
between the trapping and diffusive regimes. The decrease 
of Q!2(i) for t < t* is due to the recovery of the gaussian 
form of G s (r, t) in the trapping regime where molecules 
undergo a nearly oscillatory motion in the cages where 
they are accommodated ( fig. || ). t* provides an estimate 
of the trapping time and in fact is comparable to the a- 
relaxation time. For t « 0.2 — 0.3, corresponding to the 



second maximum of the velocity correlation function ( 
figJJ ), a small step is observed. The same feature of the 
non-gaussian parameter was observed in a MD work on 
a dense bidimensional liquid |4q] . 

The deviations of the Van Hove function by the gaus- 
sian limit at short and long r values which were observed 
for t Ri f shows that the sample has fractions of both 
trapped and highly mobile atoms 28 3^] . It has been 
suggested that the dynamics of the latter is conveniently 
described by hopping processes pgj] . To study the pos- 
sible deviations of the Van Hove function in the present 
molecular system, we considered the ratio p4|: 



N s (r) 



G a (r,t*) - Gf(r,T) 
G s (r,t*) 



(12) 



Fig. H shows the ratio N s (r) for different temperatures. 
It exhibits increasing positive deviations at both short 
and large r values on cooling. This supports the con- 
clusion that the dynamical heterogeneities evidenced in 
atomic two-phase systems are also present in molecular 
one-phase systems [plpij. 

Non-vanishing iVj^rjvalues could be anticipated by 
noting that the self-parts of the intermediate scattering 
function and the Van Hove function are related to each 
other by: 



F s (k,t) 



G s (r, t) 47rr 



sin kr 
kr 



-dr 



(13) 



According to eqjl^, the main contributions to F s (k max , t) 
come roughly from r < r* with r* w 2n/k max fa 1. If the 
discrepancies between G s (r,t*) and Gf(r, t*) for r < 1 
were small, F s (k m ax,t) « exp(— Dk^^t) for t sa t*. In- 
stead, F s (kmax j t) is found to decay as a stretched expo- 
nential ( see fig.0 and related discussion ). 

3.0 

2.0 
1.0 
0.0 
1.7 
1.2 



0.7 




2500 3000 3500 4000 
t 

1.6 - 





400 500 600 700 800 
t 



0.2 



(d) 




i ! 


T=0.520 



1100 



1200 1300 
t 



1400 



FIG. 10. Squared displacements R of selected molecules at 
T = 0.5 and T = 0.520. Note the distribution of the time 
needed to complete the jumps. 

Since the positive tail of N(r) at large r starts at r sa 1 
i.e. about the molecular size, it is tempting to ascribe it 
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to jump motion. Indeed, by inspecting the particle tra- 
jectories jumps are found. Examples are shown in fig. [l^. 
It was noted that the jump duration exhibits a distribu- 
tion ( roughly between 20 — 500 time units ), suggesting 
that different degrees of cooperativity are involved [p5[ . 
The translational jumps which are detected are relatively 
slower than the rotational jumps (flips of about 180° ) 
which take about 6 — 7 time units with little or no dis- 
tribution [ [38|]39| . Furthermore, translational jumps are 
found to be much rarer than rotational ones |3q |. In fact, 
differently from the translational case, the large amount 
of molecule flips manifests in a secondary peak of the 
proper Van Hove function p^ , p9| , pU| . 

The above remarks point to the conclusion that trans- 
lational jumps are not obviously related to rotational 
ones. Their different character will be more clearly ev- 
idenced by studying and comparing the distributions of 
the waiting-times, i.e. the lapse of time between succes- 
sive jumps ( see below and ref. @). 
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FIG. 11. Long-time (top ) and short-time ( bottom ) be- 
havior of the waiting-time distribution ij){t) at different tem- 
peratures. The solid lines superimposed to the curves is eq.|l4] 
with £ = 0.49, r = 2550 ( T = 0.5 ), £ = 0.63, r = 420 ( 
T = 0.549 ) and £ = 1, r = 49 ( T = 0.632 ). 

The evaluation of the waiting-time distribution ip{t) 
relies on the explicit definition of a jump event. In the 
present study a molecule jumps at time t if the dis- 
placement between t and t + At* (At* = 24) exceeds 
\J AR* — oaa/ 1 } = 0.5. To avoid multiple counting, 
the molecule which jumped at time t is forgotten up to 
time t + At* . Possible spurious countings due to fast rat- 
tling motion are minimized by averaging each displace- 
ment with the previous and the next one. These are 
typically spaced by 6 — 8 time units, depending on the 
temperature p^ |. It must be noted that after a time 
At* the average diffusive displacement is smaller than 
AR*, e.g. AR dlff = 6DAt* = 0.017 at T = 0.5 and 
AR di ff = 0.032 at T = 0.520. We validated the jump 
search procedure by inspecting several single-molecule 
trajectories. In particular, it was checked that molecular 
vibrations in local cages are not misinterpreted as jumps 
contributing to ip(t) at short-times. Some attempts to in- 
vestigate how ijj(t) depends on the above definition were 
made. Since they require wide statistics, the issue will 
not be touched in the present paper. It will be discussed 
in a more detailed way elsewhere. 
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FIG. 12. Comparison between the probability that the 
waiting time between consecutive jumps is larger than t, P(t), 
and F 3 {kmax,t) at T = 0.5 and 0.549. 

Fig.|ll| shows the waiting-time distribution ipit") at dif- 
ferent temperatures. At high temperature ip(t) is expo- 
nential. On cooling, the exponential decay is replaced at 
short times by a slowly-decaying regime. We fitted the 
overall decay by the function: 

i'(t) = [r(£)r«] _1 t«- 1 e- t / r < £ < 1 (14) 

The choice is motivated by the remark that in glassy 
systems rearrangements are rare events due to the con- 
straints hampering the structural relaxation. It is be- 
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lieved that intermittent behaviour in particle motion de- 
velops on cooling [ ^3||47| -|49[| . A signature of intermittence 
is the power law decay of ip(t) ( for 2D liquids see |37]] 
) and related quantities such as the first-passage time 
distribution |33| . The exponent £ of cq.[l4] has a simple 
interpretation. If a dot on the time axis marks each re- 
laxation event ( a jump ), the fractal dimension of the set 
of dots is £. For £ < 1, it is found ip(t) oc t^" 1 at short 
times p7| , p0| . Modelling the long-time decay of ip(t) is 
less obvious. As a first guess, if the distribution of events 
becomes nearly uniform ( £ = 1) beyond a time scale r 
and homogeneous across the sample, recovers the 
exponential form. 

The best fits at T = 0.5, 0.549 and 0.632 are shown 
in figjll]. The increase of temperature results in a weak 
increase of the exponent £ and a more marked decrease 
of t. It was found that eq.|lj fits ip(t) also by setting 
VAR* and At* in the ranges 0.4 — 0.7 and 24 — 48, re- 
spectively. This suggests that the character of the decay 
of ip(t) is not strongly dependent on how a jump is de- 
fined. Nonetheless, th e bes t fit values of the parameters 
£ and r depend on V AR* and At* differently. For ex- 
ample, at T = 0.5, if V 'AR* changes from 0.5 to 0.6 with 
At* = 24, £ does not change within the errors whereas 
t increases of a factor of about 2.4. r also increases by 
decreasing the time At* allowed to perform the displace- 
ment VAR*. The increase of t is understood by not- 
ing that increasing the threshold V AR* or, alternatively 
decreasing At* , reduces the number of sudden displace- 
ments which comply with the definition of "jump" and 
then increases the waiting time before a new jump oc- 
curs. Interestingly, at T = 0.5 ib(t) exhibits small but 
reproducible deviations from eq.tt4. They suggest that 
the long-time decay is faster than the exponential one. 
If the exponential decay is replaced by a gaussian one, 
the fit improves quite a lot and the £ exponent changes 
from 0.49 to 0.45. Even if this refinement is rather sug- 
gestive, we prefer to consider it as ad-hoc at the present 
level of understanding. 




Further insight on ip(t) is offered by the study of the in- 
terplay between density fluctuations and the occurrence 
of jumps. This is of interest since, according to the ex- 
tended version of MCT the dynamic transition occur- 
ring at T c from an ergodic to a non-ergodic state is actu- 
ally smeared by hopping processes [|25f Fig.[l2] compares 
at low temperatures the long-time tail of F s (k max ,t) 
and the probability that no jumps occurred before t, 
P(t) = J t °°^( x )dx. It is seen that the density self- 
correlations vanish at times when P(t) is still meaningful 
( when F s (k max ,t) ~ 0.1 P(t) ~ 0.5 at T = 0.5). There 
is also some evidence that the fraction of waiting times 
longer than the time needed to make F s (k max , t) vanish- 
ing ( e.g. F s (k max , t) < 0.1 ) increases on cooling. Strong 
analogies with the above findings are found in a recent 
MD work on the viscous silica melt |[lj . It was reported 
that the probability that a bond between a silicon and 
an oxygen atom which was present at time zero is still 
present at time t vanishes much later than F s even at 
high temperature. 

Then, at lower temperatures two restructuring regimes 
are identified depending on the length of the waiting time 
with respect to t q . This is shown in fig. |l^. For waiting 
times shorter than r Q molecules jump in a nearly frozen 
environment. Escape from the cage becomes more diffi- 
cult by lowering T and, expectedly, with a distribution of 
rates leading to the non-exponential decay of ip(t). For 
waiting times longer than t q a larger restructuring of 
the surrounding environment takes place. This averages 
the dynamical heterogeneities and leads to the long-time 
decay of ip{t). 

It is tempting to note that in the time window where 
ip(t) exhibits the power law decay MCT predicts that 
F s (k,t) itself decays as the von Schweidler power law 

S : 

F s (k,t) = fc(k)-h s (k)(-) b + ... (15) 

T 

where (fc), h B (k) and b are constants. F s (k max , t) must 
be compared to P(t) as A — Bt^ + . . . in the region of 
interest. A and B are constants. If one estimates b via 
the (3 parameter of the stretched exponential fit of the 
long-time decay of F s (k max ,t), it is found b = (3 = 0.68 
at T = 0.5. This must be compared to £ = 0.45. To make 
clearer to what extent the two time fractals are related 
to each other the analysis should be refined fiij| . This 
is beyond the purposes of the present paper. However, 
we note that for times shorter than about At* = 24 ip(t) 
and then P(t) cannot be defined since most jumps are 
not completed. 

C. Collective dynamics: the shear viscosity 

The shear viscosity 77 was evaluated by using the Ein- 
stein relation |5^]: 
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V = 



2VkT t^+oo t 



lim ~{AA(t) 2 



where: 



AA(t) = V / V aP {t')dt' 



(16) 



(17) 



VaB is one off-diagonal component of the pressure tensor 
Q,^2|j53) ( in practice AA(t) is the average over the three 
possible choices a(3 = xy,xz,yz ). 

If the quantity V a pV is evaluated based on the mo- 
tion of individual atoms comprising the molecules in the 
system we have 54 



N 



N N 

KW) y = E mtViaVip + EE f> 
= 1 i=l j>i 



'■> ) J 



(rn3-r jP ) (18) 



The sums involve components ( denoted by greek letters ) 
of Vj, r, and fy which are the velocity and the position of 
the i-th atom having mass m, and the force between the 
atoms i and j ( assumed pairwise additive ) , respectively. 
Eq.|l8| is not affected by the periodic boundary conditions 
employed in MD simulations. 

An interesting alternative to the above atomic rep- 
resentation of the pressure tensor is the molecular one 
which replaces V®g with the symmetric part of the ten- 



sor 



N N 



V«p (t)V = ]T MiViaVip + E E F ^ ( R W ~ R if>) 

i—1 i—1 j~>i 

(19) 

The sums now involve components of Vj, R,j and Fy 
which are the centre-of-mass velocity and the centre-of- 
mass position of the i-th molecule having mass Mi and 
the total force between the molecules i and j, respec- 
tively. Evaluating the atomic pressure tensor via eq.^| is 
a little more efficient than the alternative one and was 
adopted in the present study. Nonetheless, we found 
that the two representations exhibit the same conver- 
gence when evaluating eqjl^ and yield identical results 
in agreement with previous studies fl55fl. 
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FIG. 14. Temperature dependence of the shear viscosity. 
The superimposed line has slope j v ~ —2.20 ± 0.03. 

The evaluation of the shear viscosity is particularly 
time-consuming being a collective property which in- 
volves many-particle correlations. In fact, if one eval- 
uates a generic iV-particle property by a run spanning a 
time interval as long as t, the relative error r is of the or- 
der of (rc/t) 1 / 2 , t c being the correlation time. The error 
is quite larger than the one of a single-particle property 
which is of the order of 2(r c /A^) 1 / 2 Then, long 

runs are needed to reach sufficient statistics. However, 
the specific case of the viscosity is less dramatic. In su- 
percooled liquids the viscosity is roughly proportional to 
the longest relaxation time of F(k, t), the coherent inter- 
mediate scattering function, which occurs at k = k max 
ifflf . In particular, for moderately supercooled liquids a 
mode-coupling theory shows that rj is given by (57) : 



'1 



kT 
60tt 2 



dt j dkV 2 {k) 



F(k,t) 



(20) 



where V(k) = k 2 dln S{k) / dk. The vertex V(k) greatly 
reduces the weight of hydrodynamic wave vectors. The 
main contributions to the above integral are due to modes 
being located around k max with a spread \k max — fc|/A ~ 
1.4, A is the half- width of the main peak of S(k), whose 
inverse being a measure of the extent of correlations in 
direct space J57[] . Then, for a sample of volume VV the 
relative error of rj is decreased with respect to the above 
estimate by an additional factor (v/V)^ 1 ^ 2 , with v = 
A~ 3 f5f| • In order to have runs longer than the relaxation 
time of F(k max , t) and keep reasonable execution times V 
must be small. On the other hand V must be larger than 
v. We chose samples accommodating N = 108 molecules 
and carried out runs as long as 40£° at least, t° being the 
time when F s {k maxi t) vanishes (i.e. when it drops below 
0.02). t° provides an estimate of the time scale to reach 
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the limit in eqJT^. At the lowest temperature, T = 0.5, 
the volume was V ~ 150 and v ~ 8. Consequently, the 
relative error of the viscosity 2(vT c /Vt) 1 / 2 is estimated 
to be about 7%. This figure was confirmed by collecting 
several runs at each temperature. We also explicitely 
tested that the viscosity of small ( N = 108 ) and large 
samples ( N — 1000 ) at T — 0.7 exhibits no significant 
difference . 

Small samples to evaluate r\ in supercooled sustems 
were also used by Thirumalai and Mountain |28| . 

In Fig. 14 the shear viscosity is shown as a function of 
the temperature. It covers a range of more than three or- 
ders of magnitude. The data are also plotted as function 
of log(T — T c ) evidencing that the viscosity, differently 
from the diffusivity, may not be described by a power 
law analogous to eq.|^ in the overall temperature range 
investigated. If the fit with a power law analogous to 
eqj7|is limited to T < 0.85, it is found 7,, = -2.20 ±0.03. 



D. The Stokes-Einstein law 



Several experimental |3j,^,[7| |9|Jllj,|l^] and numerical 
pp| , ^4 28-3(],|33| works evidenced a decoupling of the 
translational diffusion and the viscosity on approaching 
the glass transition. Typically, the decoupling occurs 
around T c ||. To date, MD investigated the issue in 
one- and two-components atomic systems. It is therefore 
of interest to examine the present molecular system from 
that respect. 

The decoupling manifests as an enhancement of the 
translational diffusion D with respect to the prediction 
of the Stokes-Einstein relation ( SE ) which reads 1 58 



(21) 



fi is a constant that depends on both the molecule ge- 
ometry and the boundary conditions. For a sphere of 
radius a, fj, equals Qira or Ana if stick or slip boundary 
conditions occur, respectively. The problem of uniaxial 
ellipsoids in the presence of stick boundaty conditions 
may be worked out analytically |S8| . Tables of fi for pro- 
late ellipsoids with slip boundary conditions were also 
reported J59| . The case of the biaxial ellipsoid with stick 
boundary condition was discussed recently by noting an 
interesting electrostatic analogy j60|. 
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FIG. 15. The temperature dependence of the ratio Dr//kT. 
SE predicts a constant value. Dashed lines are the SE predic- 
tions for prolate ellipsoids with semiaxis b = 0.46 and c = 0.69 
and stick or slip boundary conditions. A magnification of the 
low-T region is shown by the inset. 

In Fig. [ij] we plot the ratio Drj/kT as a function of 
temperature. According to SE the ratio must be con- 
stant. At higher temperatures the ratio levels off at about 
0.105 ±0.007. On cooling, there is first a mild change fol- 
lowed by a steep increase below T = 0.632 = 1.38T C . 

It is believed that the SE failure is a signature of the 
heterogeneous dynamics of supercooled liquids [|l6] [l9| . 
Alternative views are provided by frustrated lattice ga s 
models |20) and the "energy landscape" picture ||,^l]-||4| . 
Most intepretations suggest the existence of crossover 
temperatures broadly located around T c ||. From this 
respect it is tempting to note that the SE law breaks 
down at T ~ 0.632 below which intermittent behavior ev- 
idences ( see fig.fn] ). Intermittence disappears for times 
longer than r a ( see fig. |lj ) which is consistent with a 
growing dynamic heterogeneity of the liquid. 

Around T = 0.77 a plateau at 0.151 ± 0.01 is reached. 
We have compared the results to the prediction of the SE 
law for prolate ellipsoids. The diatomic molecule under 
study may be roughly sketched as a prolate ellipsoid with 
semiaxis 6 = 3/4 and c = 1/2. The corresponding ratio 
Drj/kT for stick and slip boundary conditions is equal 
to 0.091 and 0.1415, respectively. The values compare 
well to the high- and low-T plateau which are observed 
in fig]!?]. By setting b = 0.69 and c = 3/26 the agreement 
is improved with Drj/kT = 0.098 and 0.154 for stick and 
slip boundary conditions, respectively. The above anal- 
ysis provides reasonable evidence of a precursors effect 
of the SE breakdown which manifests itself as an appar- 
ent stick-slip transition. A quite similar crossover from 
stick to slip boundary conditions has been observed on 
approaching the glassy freezing of colloidal suspensions 
|6l| . Since we are studying a one-component system, the 
effect seems to be conceptually different by the apparent 
change of the boundary conditions which has been evi- 
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denced by a recent MD study of the motion of a guest 
tracer in a liquid host 0] . 

It is worthwhile to mention that the SE law evi- 
dences hnc modifications of the supercooled liquid behav- 
ior which are hardly found by inspecting the diffusivity. 
The latter exhibits a single power-law regime over all the 
temperature range we studied. 



product Dri/T increases and the Stokes- Einstein relation 
is not obeyed. The breakdown occurs near to the tem- 
perature where the intermittency is evidenced by ip(t). 
Interestingly, a precursor effect of the breakdown is ob- 
served which manifests as an apparent stick-slip transi- 
tion. A crossover from stick to slip boundary conditions 
has been observed on approaching the glassy freezing of 
colloidal suspensions . 



IV. CONCLUSIONS 

The transport and the relaxation properties of a molec- 
ular supercooled liquid on the isobar P = 1.5 has been 
studied by molecular dynamics. The molecule is a rigid 
A-B system. On cooling, the diffusivity decrease is fit- 
ted over four orders of magnitude by the power law D oc 
(T - T C ) 1D with 7i3 = 1.93±0.02 and T c = 0.458±0.002. 
The divergence of the primary relaxation time r a drawn 
by the intermediate scattering function F s (k rnaxi t) tracks 
the diffusivity at the lowest temperatures according to 
r a oc (^ OI fl) _1 . The result confirms findings on hard 
sphere systems [EEJ. It disagrees with other studies on 
the same model system under study here which noticed 
that a power-law fit of D and r Q yields the same T c value 
but la > lD with N = 500, P = l§. 

At the lowest temperatures fractions of highly mo- 
bile and trapped molecules are evidenced, then extend- 
ing previous results on supercooled atomic mixtures to 
one-component molecular liquids ]28| , |3~4| ] . Translational 
jumps are evidenced. The duration of the jumps ex- 
hibits a distribution. The distribution of the waiting- 
times before a new jump takes place, ip(t), is exponen- 
tial at higher temperatures. At lower temperatures two 
regimes are evidenced: at short times ^ ^ with 
< £, < 1 whereas at long times the decay is faster than 
exponential. The crossover between the two regimes oc- 
curs around r Q . Noticeably in this time window MCT 
predicts that F s (k max , t) also decays as a power law. The 
interplay between the two time fractals is anticipated and 
a preliminary analysis has been attempted. The fractal 
distribution of the waiting time is ascribed to the inter- 
mittent behavior which is expected to develop in glassy 

"-HJH 



systems |4' 
exceeds r„ 



If the waiting time before a new jump 
the environment surrounding each molecule 
largely restructures. The subsequent average process re- 
sults in a faster decay of i/j(t). The probability that no 
jump occurred before time t is found to vanish slower 
than F s (k max , t) in close analogy with the case of viscous 
silica melts 

The shear viscosity has been studied over more than 
three orders of magnitude. The data are fitted by the 
power law 77 oc (T — Tc) 7 " with 7,, = -2.20 ± 0.03 at the 
lowest temperatures. The validity of the Stokes-Einstein 
relation has been examined. Previous MD studies were 
limited to atomic one- and two- components systems. At 
higher temperatures SE fits well the data if stick bound- 
ary conditions are assumed. At lower temperatures the 
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APPENDIX A: RATTLE-NPT ALGORITHM 

The well-known Nose-Andersen algorithm ensures 
NPT equilibrium, i.e. equilibrium at constant number of 
particles, pressure P and temperature T GJ|. It defines 
one extra degree of freedom s(t), describing the thermal 
bath and allows the change of the volume V(t). Since 
the algorithm was originally derived for atoms, the ex- 
tension to molecules is of interest. To this aim, the so- 
called constraint methods are of help, particularly the 
RATTLE algorithm KJ|. It evaluates the dynamics of 
polyatomic molecules by defining proper forces constrain- 
ing the relative motion of the atoms belonging to the 
same molecule. Below, the combined RATTLE-NPT al- 
gorithm is described. 

First, the vector R(t + St) of all atoms positions, 
s(t + St) and the vector V(t + 5t) of all atoms veloc- 
ities and the their first time derivatives at mid-step 
V(t + St/2,s(t + St/2) , V(t + St/2) are evaluated ac- 
cording to the RATTLE algorithm . The forces f;(£ + St) 
(i refers to the i-th atom) are expressed in terms of the 
configuration H(t+5t). The other steps are the following: 

1. Guess of the first derivatives of interest at time t + 
St, y°(t + 5t) (y = v,s,V): 



y 9 (t + 5t) = 2y(t)~y(t-St) 

The above guess, which is correct to second order, 
is not unique and alternatives are possible. 

2. Calculation of temperature T(t+St) by using v 9 (t). 

3. Calculation of s(t + St). Replacing the lagrangian 
equation for s(t + St) into the velocity Verlet equa- 
tions Q yields: 
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s(t + St) = A + ^-St + ^St 2 + 0(St 3 ), 

with ,4 = s{t+\5t) + \5tf T[t+5 g- To s{t+5t). Q, T 
and / are the mass of the thermal piston, the ther- 
mal bath temperature and the degrees of freedom 
( / = r;N at — 3 in the present case ), respectively. 

4. Calculation of pressure P(t + St). 

5. Evaluation of V(t + St) according to the velocity 
Verlet algorithm: 

V(t + St) = 
y(t + \5t) + lity(P(f + St) - Pp) 
1 - ±6U(t + 5t)/s(t + 8t) 

where W and P(t) are the mass of the piston setting 
the pressure at P and the instantaneous pressure, 
respectively. 

6. Evaluation of the velocities at time t + St, Vi(t + St) 
in the presence of the intermolccular forces, the 
thermal bath forces and the forces due to the me- 
chanical piston. According to the velocity Verlet 
algorithm: 

Vi(t + St) = (1 + \Sts{t + 5t)/s(t + St))- 1 x 
\v i (t+Ut) + ytf i (t + 5t) + 



hst 




R c 



where R cm ,i is the center of mass of the molecule 
where the i-th atom is located. 

7. Adjustement of the velocity Vj. The relative ve- 
locities of atoms which are mutually bonded are 
considered and their component along the bond is 
made to vanish, according to the RATTLE algo- 
rithm. 

Finally we note that iterating steps from 2 through 7 
increases the accuracy considerably. The present algo- 
rithm needs to store the arrays V(t — St), s(t — St) and 
V(t — St) to guess the velocity as v 9 (t + St) in the first 
step. 
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